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Abstract The article considers symmetric general linear methods, a class of nu¬ 
merical time integration methods which, like symmetric Runge-Kutta methods, 
are applicable to general time-reversible differential equations, not just those de¬ 
rived from separable second-order problems. A definition of time-reversal symme¬ 
try is formulated for general linear methods, and criteria are found for the methods 
to be free of linear parasitism. It is shown that symmetric parasitism-free methods 
cannot be explicit, but a method of order 4 is constructed with only one implicit 
stage. Several characterizations of symmetry are given, and connections are made 
with G-symplecticity. Symmetric methods are shown to be of even order, a suitable 
symmetric starting method is constructed and shown to be essentially unique. The 
underlying one—step method is shown to be time-symmetric. Several symmetric 
methods of order 4 are constructed and implemented on test problems. The meth¬ 
ods are efficient when compared with Runge-Kutta methods of the same order, 
and invariants of the motion are well—approximated over long time intervals. 

Keywords time-symmetric general linear methods • G-symplectic methods • 
multivalue methods ■ conservative methods 
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1 Introduction 

Symmetric general linear methods are a class of multistage multivalue methods 
with time-reversal symmetry. As we demonstrate, such methods can efficiently in- 
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tegrate the solutions of differential equations which are themselves time-reversible, 
in such a way that invariants of the motion are preserved over long time inter¬ 
vals. The main aim of this paper is to characterize, construct and test high-order 
symmetric general linear methods with minimal implicitness and zero parasitic 
growth-parameters. 

Under mild conditions, the flow associated with a general ordinary differential 
equation satisfies the basic time-reversal symmetry E-xEx = I. A Runge-Kutta 
method is symmetric if it satisfies the analogous property. 


M_hMh = I, 


( 1 ) 


where M/j is the map generated by a single step of the method. As shown in |24| . 
m, o is also sufficient for a Runge-Kutta method to inherit the stronger sym¬ 
metry of p-reversibility, when a differential equation has this property. Symmetry 
implies even order and leads to simplifications in the order theory for such meth¬ 
ods, [19]. Practically, symmetric Runge-Kutta methods are shown to perform well 
for such problems over long time intervals in the book of Hairer, Lubich & Wan¬ 
ner m- However, every irreducible stage of a symmetric Runge-Kutta method 
is necessarily implicit, |22| . |28|.|10l V.2]. The most efficient such methods are 
DIRKs, formed by compositions of the implicit midpoint method, ED, m, m, 
|18| . For separable problems originating from a system of second order differential 
equations, the symplectic Euler and Runge-Kutta-Nystrom methods have been 
generalized to obtain higher order partitioned Runge-Kutta methods m , some 
of which are explicit. The most popular low order method for separable problems 
is the explicit Stormer-Verlet method ED, which may be viewed as a partitioned 
Runge-Kutta method, a partitioned linear multistep method, or a non-standard 
implementation of the leapfrog method. 

The properties of standard linear multistep methods and one-leg methods were 
investigated by Eirola & Sanz-Serna |6] , who showed that symmetry is equivalent 
to G-symplecticity in this case. The properties of symmetric multistep methods 
were further investigated in |3]. However, Dahlquist |5] had already shown that 
the parasitic roots of such methods have non-zero growth-parameters. Hence, 
symmetric linear multistep and one-leg methods are weakly unstable. 

An important class of systems with time-reversal symmetry are of the form 


fy 

dx^ 


= /(y). 


familiar from many examples in Mechanics and other branches of Physics. The clas¬ 
sical Stormer-Cowell linear multistep methods , popular with Astronomers, 

exploit the special structure of such systems by directly approximating the sec¬ 
ond derivative. However, only lowest order method (Stormer-Verlet) is symmetric. 
The articles |5] and m made early studies of the stability properties of second 
order multistep methods. New symmetric high order second order methods were 
designed and successfully tested in [20]. Hairer & Lubich 0, ED used backward 
error analysis techniques to show that the underlying one-step method is a sym¬ 
metric approximation of the true solution, and that parasitic components remain 
under control for long times. 
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As a model for general linear methods, consider a fc-step linear multistep 
method in one—step form. Here, M/j may be interpreted as the map 

[Vn, • • • , Vn+k-lj hfn, • ■ ■ , 

n £ L. (I) 

' [l/n+l5 ■ ■ ■ 5 yn-\-k^ • • • 5 

Under the change of variable m = n + k, represents the mapping 


[ym+fc —l! ■ ■ ■ 1 Vm, /t/m-j-Zc—15 ■ • • 5 kfm] 

[Vm+k^ • • • 5 Vm+l: ■ 5 


m € Z. 


To equate this mapping with a coordinate transform L is needed, which re¬ 
verses the order of both the y and hf entries. Furthermore, L must multiply the 
hf terms by —1. (Both these actions of L are directly related to time-reversal.) 
Then, the following modification of m holds: 

L'M_hL'Mfi = I, such that = I. (3) 


As shown in Section|3l identity (IS]) characterizes a symmetric general linear method. 
In Section [S] it is shown that ((S]) implies that ([!]) is formally satisfied by the cor¬ 
responding underlying one-step method. These results are essentially similar to 
those of nni XIV.4.2], though our assumptions differ in detail. 

Three further characterizations of symmetry are obtained in the paper: 

(i) In Section |3l an algebraic condition (EOl) in terms of the method coefficient 
matrices (A, U, B, U), the matrix L and a stage permutation matrix P, which also 
satisfies = I\ see also [M XIV.4.2]. This condition, together with the canonical 
form identified later in Section [S] is the most useful in method construction. 

(ii) In Section m an AA-stability condition: LM{—PZP)LM{Z) = I for all suf¬ 
ficiently small diagonal Z, where M{Z) is the non-autonomous linear stability 
matrix, cf. [T] . This condition helps to show linear stability on a subinterval of the 
imaginary axis. 

(iii) Also in Section 31 a characterization in terms of the matrix transfer function, 
generalizing the one-leg condition of [6], (cr/p)(C) = —{a/p){(~^). This condition 
has potential application to long-time nonlinear stability theory, and also helps in 
the construction of methods that are both symmetric and G-symplectic. 

Parasitism is a potential disadvantage for any non-trivial symmetric general 
linear methods. However, in Section 31 we find necessary and sufficient conditions 
on the coefficient matrices of the method for the linear stability matrix M{zl) to 
have sublinear growth in parasitic directions. (In the terminology of |10l XIV.5.2], 
this is equivalent to all parasitic roots having zero growth-parameters.) These 
coefficient conditions play a critical role in the construction of practical methods 
in Section |6l They are also used to show that there are no explicit symmetric 
parasitism-free methods. 

In Section [5l it is shown that a symmetric general linear method is always 
of even order. Central to this result is a constructive proof of the existence and 
uniqueness of a symmetry-respecting starting method S/j satisfying 




( 4 ) 
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with respect to which M/j is of maximal order. Related ideas are used to show the 
existence of a formal starting method, underlying one-step method pair <Pfi) 
such that 

= ^h^h- (5) 

Example symmetric methods of order 4 are constructed in Section |6l These 
methods have diagonally implicit stage matrices, and some are also G-symplectic. 
One method has only one implicit stage, and is therefore theoretically more ef¬ 
ficient than a symmetric DIRK of the same order. The simulations of Section 0 
show that symmetric general linear methods approximately conserve the Hamil¬ 
tonian of several low-dimensional symmetric problems over long time intervals in 
a similar way to symmetric Runge-Kutta methods. Furthermore, there are 4th 
order symmetric general linear methods with fewer implicit stages than is possible 
in the Runge-Kutta case. This leads to some efficiency savings over long-times. 


f X ^ X, and yo ^ X, let y = yyg denote the solution of the 

( 6 ) 


2 General linear methods 

For X = 

autonomous initial value problem, 

y'{x) = f{y{x)), y(0) = yo- 

For a; G M, denote the flow for hy Ex ■ X —> X, so that 

y(x) = Exyo- 

For all ODEs, the evolution operator satisfies the group properties, 

Eq = /; Ex-^Ex2 = 3 ;i, X2 G K; Ex = E—x, a; G K (7) 

We refer to a general linear method (H, U, B, V), where 


A U 
B V 


(8) 


forms a partitioned (s-|-r) x (s-)-r) complex-valued matrix or tableau. For practical 
methods, the coefficients are real, but for some theoretical purposes the complex 
case is also treated. 

For time-step h and n G N, the new values G X^ are found from G X'^ 

via the formulae 


Y = h{A® I)E +{U® 
= /i(R (g) 7)F + (K (g) 


( 9 ) 

( 10 ) 


defined using temporary Y,E £ X®. The subvectors in F (the stage derivatives) are 
related to the subvectors in Y (the stages) by Fi = fiYi), i = 1,2,... ,s. Usually, 
where no ambiguity is possible, the Kronecker products in m and m will be 
omitted and we write 

Y = hAF + Uy^'^~^\ 


= hBF + Vy 


[n-l] 


In this paper, the method is always assumed to satisfy the conditions below. 
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Definition 1 A general linear method {A,U, B,V) is 
(Pi) Preconsistent, if {1, u, w”) is an eigentriple ofV, sueh that 

Vu = u, = w”, w^u = 1; 

m Consistent, if it is preconsistent, Uu = 1, and there exists non-zero v € C’~ such 
that B1 + Vv = u + v; 

(P:) Zero-stable, if l|l^"ll < oo. 

To approximate the solution of (ED with initial data yo £ X, we generate 

= Syyo, 


using a practical starting method Sy : X —> , where the tablean 

'A l' 

B u ’ 


( 11 ) 


has dimensions (s + r) x (s + 1). Similarly, a practical finishing method, : X’' 
X, is reqnired. It is assnmed that Ff^o = Ix- 


3 Symmetric methods 

We define symmetry in the context of the nonlinear map generated by the method. 
Other characterizations of symmetry are considered, with a view to identifying or 
constructing symmetric methods. 


3.1 The method as a nonlinear map 


For f : X ^ X and time-step h, the method maps an input vector y G X^ to an 
output vector M/ji/. Define the nonlinear map M/j : X'" —> X^ by 


Y = hAF + Uy, 
Myy = hBF + Vy. 


( 12 ) 

(13) 


(It will be assumed that / and h are such that (l9|) has a solution, and that a 
suitable selection principle chooses a unique Y when multiple solutions exist.) 
Equivalent maps: The map M/j is not changed if a different ordering is chosen 
for the subvectors of Y ; that is, M/j is also generated by the method defined by 
the tableau 


P-^AP P-^U' 
BP V 


(14) 


where P is a permutation matrix. 

If r G ([^rxr non-singular, then T~^JAhT is equivalent to M/j. The identity 


{T-^MhT){T-\) = T-^Mhy, 


shows that T only changes the coordinate basis. A tableau for T is 


A UT 
p-^B T~^VT 


( 15 ) 
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3.2 Symmetry of the map 

We say that the map is symmetric if the process of calculating from can 
be reversed by using an equivalent map with the sign of h reversed; i.e. 

Mh = LMzIl, (16) 

for some nonsingular matrix L G such that = I. Physically, the involution 

L corresponds to a linear change of coordinates for y to take account of the change 
in time direction. Algebraically, the condition = I is required to ensure that we 
recover after two iterations of (nnD. This definition is similar to that stated in 
[?, XIV]. 

The inverse map: From m and m, we deduce that the inverse map ^ 
satisfies 


Y = hAF + UMZ^y, 
y = hBF + VMZ^y. 

Solving these equations for yields 

Y =-h{UV~^B - A)F+ UV~^y, (17) 

MZ^y =-hV~^BF+ V~^y. (18) 

3.3 Symmetry of the method 

We say that the method is symmetric if 

'AT PAP - UV~^B PU - ULV 
BP - VLB L - VLV 

More specifically, we say that method (A, U, B, V) is (L, P)-symmetric if (1191) 
holds. 

Proposition 2 Suppose that JAh is the map associated with a symmetric method. 
Then, JAh is symmetric. 

Proof A rearrangement of definition dUD yields 


= 0, = 7, = 7, (19) 


'A U' 


'P(UV~^B - A)P PUV~^L' 

B 


LV~^BP LV~^L 


( 20 ) 


Here, the left-hand side of dSOl) is a tableau for JA^. Taking note of (HI, (nsi), 
(ED and (ED, the right-hand side of (EOD is one possible tableau for 7,M_^7,. 

Remark: The tableau on the right-hand side of [20] is also known as an adjoint 
tableau for the method (A, U, B, V). The conditions = I and P^ = I ensure 
that the original tableau is recovered after 2 iterations of (pop. The coefficient 
conditions in (1201) are similar to those given in [101 XIV], except that L and P are 
not involutions there. 
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Fig. 1: Relationships between various mappings 



Fig. 2: The role of the underlying one-step method 


3.4 Symmetry of the starting method 

In order to ensure that MJjS'/i = it is required that the starting method 

satishes 

S_h = LSh. ( 21 ) 

Considering the tableau (HH, this is equivalent to the coefficient conditions 

A = -PAP B = -LBP, (22) 

for some permutation matrix P £ such that P^ = 7. 

The diagram in Figure [1] shows the relationship between various quantities and 
mappings which have arisen in this discussion. In addition to 'Mh^ we introduce a 
further mapping defined as T = in Q. 

In Figure[2l the role of the underlying one-step pair (S^, 45^), discussed in the 
Introduction, is also included. 


3.5 Canonical form based on F-diagonalization 

Given a stable consistent general linear method {A,U, B,V), which is {L,P)- 
symmetric, we explore a canonical form of the method based on a diagonal form of 
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V. The approach is to successively transform {A, U, B,V) to an equivalent method 
{A,UT,T~^B,T~^VT) and then to regard this as the base method. This leads 
to a specific form for the coefficient matrices of the method which can then be 
back-transformed to a convenient format for practical considerations, such as a re¬ 
quirement that U, B, V should be real matrices. As transformations to a canonical 
form take place, L is also transformed. 

Methods in canonical form are convenient to analyze in terms of order of ac¬ 
curacy and the possible presence of parasitic growth factors. 

Since V is similar to V~^ and each is power-bounded, T exists such that T~^VT 
is diagonal with diagonal elements made up from points on the unit circle. We will 
see how to carry out this diagonalization process in such a way that, when the 
corresponding transformation has also been applied to U and B, these matrices 
have a specihc structure. Because of the original real form of V , the diagonal 
elements are real or come in conjugate pairs. Hence, in the canonical form, 

H = Ho © T-i © Vi © V2 © ■ • ■ © 

where 

Ho = diag(l, 

H_i =diag(-l, 

H = cliag(^j, (^ 2 ,..., Cii • • • 5 f = 1, 2,..., n. 

The number of diagonal elements in these blocks are respectively mo, m_i and 2mi. 
Because of consistency of the method mo > 1, but it is possible that m_i = 0, 
indicating that this block is missing. It is assumed that m^ > 1, although it is 
possible that n = 0 indicating that the hnal n blocks in H do not exist. 

To carry out the diagonalization process, dehne transforming matrices T and 
T~^ of the forms 

5-0 

^-i 

Si 


where, the various submatrices are blocks of eigenvectors; that is 

HTo=To, HT_i = -T_i, VT, = QT„ VTi = i = l,2,...,n, 

SoV = So, S_iH = -S_i, S,V, = QS„ S,V = C,Si, i = l,2,...,n, 
with To, r_i, So, S_i real. Similarly, transformed B and U matrices have the form 

U= [Uo U-i UiUi--- UnUn]. 



r= [To r_i Ti Ti ••• Tn T„] , T-1 = 


B = 


(23) 
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In the canonical form, V ^ = V and we recall that L ^ = L, so that (EOl) 
becomes 

LV = VL, 

and it follows that, in a block representation of L compatible with the block 
structure of V, the off-diagonal blocks are zero. Hence we can write 

H = Lq 0 L — i 0 Li 0 Z/2 0 • ■ ■ 0 Ln- 


We now consider the structure of the diagonal blocks in L. In the case of Lq and 
L-i, the idempotent property implies that these matrices are similar to diagonal 
matrices of the form 7 0 (—7) where the dimensions of the +7 blocks and the —7 
blocks are not necessarily the same. Hence, by imposing additional transformations 
on the method if necessary, we can assume this diagonal form for Vo and V-i. For 
the blocks Vi, i = 1,2,... ,n, the equation LiVi = ViLi implies that there exists a 
non-singular x matrix Ki such that 


L^ = 


0 K, 
K-^ 0 


(24) 


The choice of the non-singular matrix Ki is arbitrary. To see why this is the case, 
apply the transformation 


■7 0 ■ 

T/ 

■7 0 

0 K, 

Vi 

.0K-\ 


which leaves V) unchanged. The transformation (I25II applied to Li gives 


(25) 


■7 0 ■ 


0 K' 


■7 0 


0 7' 

0 7^, 


0 


.0K-\ 


7 0 


so that Ki has been replaced by 7. We will take the canonical form of Li to be 
(1241) with Ki = 7. 

Using the new basis, with Ki = I gives 

7/ I —>UT and Hi—^ T~^B, 


where, as indicated above, we now use U and B for the transformed matrices. A 
rearrangement of the symmetry conditions dMI now yields 

B = VLBP, U = PULV. (26) 

Using the canonical forms of V and L, and taking j >1, (12611 implies that 

[Uj, = [PUj, PUj] 


'Bj' 


0 

Cji' 

'BjP' 

[Bj\ 


.cy 

0 

[BjP\ 


0 Cj7 
C,7 0 


for the submatrices Bj G Uj € . This simplifies to 

Bj=CjB,P, U,=CjPUj. (27) 

If P represents the permutation tt, then the components of Bj and Uj satisfy 

^ji ^ij 1 < Z < S. (^8) 
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3.6 Formulation in real form 


Having constructed a method in canonical form, it is desirable to transform it back 
to a formulation in which B,. U and V have only real elements. Consider complex 
blocks Bi,Bi and Ui,Ui in ([23ll. corresponding to Vi = diag((^/, ^/). We will show 
how it is possible to construct Ti so that 


T~ 


B, 

B, 


[u,u^]T, 71-Vr 


are each real. The suggested choice of Ti and 71 ^ are 


7 if 

, T~^ = i 

7 f 

I -il 

’ 1 2 

-il il 


leading to transformed blocks 


'Bf 


Re Si 



Im Bi 


[[/, U,]T= [2ReU, -2Imt/,] , 


7i-idiag(C,7,Ci/)r, 


ReCil - ImQI 
-ImCiJ ReC7 


4 Stability 

4.1 Linear stability 

Definition 3 For a method {A, U, B, V) and Z = diag(j 2 i, . .., Zs) € such that 

I — AZ is non-singular, the linear stability function is given by 

M{Z) :=V + BZ{I - AZ)~^U. (29) 

Theorem 4 Method {A, U, B, V) is symmetric, if and only if there exists a permuta¬ 
tion matrix P with = I such that for all diagonal Z £ with ||yl||||Z|| < 1, 

LM{-PZP)LM{Z) = L (30) 

Proof (only if) Assume first that A is B-irreducible, see |15| . Choose t/[°l £ C® and 
diagonal Z £ C®^® such that ||A||||Z|| < 1. Let Y be the unique solution of 

{I-AZ)Y = 

For almost all Z, S irreducibility implies Yi ^ Yj implies i ^ j, 1 < «, i < s. Using 
interpolation, we may construct continuous / : C —> C such that 

fiVi) = ZiiYi, 1 <i< s. 

Let h = I and set := BP + Vy^^^ Then, 

yW =M(Z)y[°l. 
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Also, by (EOl), = LM{-PZP)Ly^^\ Thus, 

y[°l = LM{-PZP)Ly^^'^ = LM{-PZP)LM{Z)y^°\ 

Hence, (EOl) holds for almost all diagonal Z, ||Z||||A|| < 1, when A is S'-irreducible. 
The general case follows from the continuity of LM(—PZP)LM(Z) and the density 
of Z and S'-irreducible A. 

(if) Sketch Let Z = t diag(x) for non-zero a; € C'*, and expand (1301) in powers of e. 
The identities in the 4 quadrants of (EOl) follow from the terms of 0(1), 0(e) and 
0 (e2). 

Remark: This result, which generalizes a linear multistep theorem of [6], is in 
the spirit of the AA'-stability characterization of algebraic stability [1]. In the 
Runge-Kutta case, L = 1, and dsni) generalizes the known necessary condition for 
symmetry: M{—zI)M{zI) = I |10l V.6]. 

Lemma 1 Suppose that the method {A, U, B, V) is (L, P)-symmetric and has real 
coefficients. Suppose also that i? G is a diagonal matrix such that R = PRP. 

Then, 

C € a{M{iR)) C”^ e cr(M(iR)). 

Proof The symmetry of M and the assumption on R imply that 

LM{iR)LM{-iPRP) = I. 

Hence, M{iR) is of full rank and possesses an inverse. In particular, ( ^ 0 and 
G cr(M(i7?)“^). Since 

LM{iR)~^L = M{-iPRP) = M(-iR), 

similarity implies that G a{M{—iR)). Taking the complex conjugate, it follows 
that e o{M{\R)). 

Theorem 5 Assume that the method (A, U, B, V) is symmetric and has real coeffi¬ 
cients. Assume also that the eigenvalues of V are distinct. Then, there exists ko > 0 
such that the eigenvalues of M{iR) are distinct and unimodular for all diagonal R G 
such that R = PRP and ||i?|| < ko. In particular, the linear stability domain S 
contains an imaginary interval (—ifco, ifeo)- 

Proof The eigenvalues of V are unimodular, so C £ cr(H) implies ( — C ^ = 0. Let 
(5 > 0 be the closest distance between any two eigenvalues of V. For small diagonal 
i? G K'*, the eigenvalues of M{iR) are continuous functions of R. Thus, there exists 
fco > 0 such that ||i?|| < ko implies 

(i) no eigenvalue of M{iR) is closer than 5/2 to any other eigenvalue; 

(ii) C £ o-{M{iR)) implies |C — C ^1 < ^/4- 

Now, suppose that diagonal R G satishes ||i?|| < ko and that (/ G a{M{\R)). 
Then, Lemma [1] implies that ^ G a{M(iR)). Furthermore, conditions (i) and (ii) 
imply that (( = i^; i.e. |C| = 1- 

If i? = xl, \x\ < ||r7|| < ko, then the foregoing results imply that all eigenvalues 
of M{ixl) are unimodular. Thus, M{ixl) is power-bounded and ia: G S. 

Remark: The continuity argument used in the proof of Theorem [5] may be used 
to increase ko until M(ifco) is ill-dehned or has a multiple eigenvalue. 
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4.2 Parasitism 

For a zero-stable symmetric method, V is power-bounded and similar to V~^. 
Hence, all of the eigenvalues of V are unimodular and, at worst, semi-simple. 
Typically, however, symmetric methods will be applied to problems without overall 
growth or decay. Hence, care is needed to limit the growth of components of the 
numerical solution associated with the non-principal eigenvalues of V. Below, it 
is assumed that the method is written in the canonical coordinates of Subsection 
[331 

Definition 6 A preconsistent symmetric method is said to be parasitism—free if there 
exist C, v > 0 such that, given e > 0, 

\\{{I-eiel)M{zI){I-eiel)r\\<C(l + ue^)^, n e N, (31) 

for all z € C such that \z\ < e. 

Proposition 7 A preconsistent symmetric method is parasitism-free if and only if 

wlBUu^ = 0, (32) 

whenever C, is a non-principal eigenvalue of V, and and are respectively right 
and left eigenvectors corresponding to 

Remark: In the assumed canonical coordinates, (1321) implies the simple condition, 

{BU)u = 0, 2<i<r, (33) 

If ^ is a multiple eigenvalue of V, (13211 implies that some off-diagonal elements are 
also zero. 

Proof If (1321) holds, then there exist r eigentriples (C, U(^, w”) for V such that 
VU(^ = 

Set T = [u^j, ..., Then, T~^ = ..., and T~^VT = diag(Ci, •. ■, Cr)- 

Assuming is the principal eigenvalue, 

T-\l - eiel)M{zI){I - eiel)T = T-\l - eiel){V + zBU + 0{\z\^)){I - eWi)T 
= diag(0, C 2 , ..., Cr) + - eiel)BU{I - eiel)r + 0(\zf). (34) 

For small z € C, eigenvalue perturbation theory (Wilkinson 1965) implies that the 
eigenvalues of {I — eie\)M{zI){I—eie\) consist of a term of 0(|2:|^), corresponding 
to the principal eigenvector of V, and {(j{z)}'^^ 2 ^ where 

0 («) = 0(0) + zwlBUu^. + 0{\z\^) = 0(0) + O(0p), 2<j<r. 

Since 10(0)1 = 1, the parasitism-free condition (13111 is satisfied. 

Conversely, assume that (|32ll is not satisfied, and let small e > 0 be chosen. Set 
z = 6C/w'^BUU(^, where 5 > 0 is such that \z\ = e. Then, similarly to (13411 . 

wq((I - eiel)M(zI){I - eiel))" m^ = C”(1 + ne + 0(e^)), 
and so m cannot be satisfied. 
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Corollary 1 There are no explicit consistent symmetric parasitism-free methods. 

Proof Consistency o, with u = w = ei, implies that 

{BU)ii = elBUei = elBl = el(7 - V)v + ejei = 1. 

Combined with (f33l). the canonical form for V = diag(l, ( 2 , ■ • •, Cr)j and the first 
quadrant of (1201) . we deduce that 

r 

1 = Y^{BU)uC~^ = tr{BUV~^) = ti{UV~^B) = tr(^ + PAP). (35) 

i=l 

Hence, diag(j4) ^ 0. Thus, the method has at least one implicit stage. 

Remark: As mentioned in the Introduction, it is known that all symmetric Runge- 
Kutta methods are implicit, and that all symmetric linear multistep methods suffer 
from parasitism, (whether or not they are explicit). An example in Section [6] 
show that only one implicit stage is necessary for a general linear method to be 
symmetric and parasitism-free. 


4.3 Transfer function characterization of symmetric methods 


Definition 8 For a method (A, U, B, V), and G C such that C,I —V is nonsingular, 
the transfer function is defined by 

N{()■.= A PUiCI-V)~^B. (36) 

This function has previously been considered in [1] and m in the context of 
algebraically stable methods. We omit the proof of the following straightforward 
result. 


Lemma 2 ([T]) Given GLMs {A, U, B, V) and {A, U, B, V), with diagonalizable V 
and V, 

N(0 = N{0, CeC\{a{V)Ua{V)u{0}), 

if and only if there exists non-singular T G (prxr 


A U' 


'A UT ' 

.B T. 


p-ip T~^VT 


(37) 


Theorem 9 A method {A, U, B, V) is symmetric if and only if there exists a permu¬ 
tation matrix P such that P^ = I and 

N{C) = -PN{C^)P, CeC\Ao, (38) 

where Aq : = {C £ C : |C| = 1 or ^ = 0}. 

Remark: Identity (1381) is an L-free characterization of symmetry, which generalizes 
the {o-/p){() = —(cr/p)(C“^) condition [6] for multistep symmetry. 
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Proof (only if) Given a method {A, U, B, 1/), the method (A*, U*, B*, V*) ap¬ 
pearing on the right-hand side of (poll is the adjoint method, see m- From formula 

N*(C) = P{UV~^B - A)P + PUV~^L{CI - LV~^L)~^LV~^BP 
= P{-A + U{I -h {(V - I)~^)V~^B)P 

= -P{A + U{C^I - V)~^B)P = -PN{C^)P, (39) 

for C € C \ Z\o- For a symmetric method, (12011 implies that {A*, U*, B*, V*) = 
{A, U, B, V). Hence, (l39l> implies (l38ll . 

(if) Now assume that (I38II holds for method [A, U, B, V), and let {A*, U*, B*, V*) 
denote its adjoint for L = I. From identity (1391) . we know that 

n{Q = -pn{c^)p = n*{Q, CeC\zio. 

Applying Lemma [2] there exists non-singular T € 


'A U' 


' P{UV~^B - A)P PUV~^T' 

B V_ 


T~^V~^BP T~^V~^T 


Using a diagonal decomposition of V, as in Subsection 13.51 T may be altered if 
necessary so that = 7 on each eigensubspace of V, without affecting identity 
(|i0l). Thus, (001) holds for T = L auch that = I. 


4.4 A transfer function characterization of G-symplectic methods 


It is the purpose of symplectic, or canonical, one-step methods to preserve the 
value of [yn,yn\Q as n increases, where the symmetric bi-linear function [•,-]q is 
defined by 

[y,z]Q ■= {y,Qz}, 

Q is a symmetric N x N matrix, and (•, •} is an inner product on A = If 
[y,f{y)]Q = 0, then [y{x),y{x)]Q is an invariant of the ODE (l6|). 

For a general linear method (l8|), it is necessary to work in the higher dimen¬ 
sional space and we consider the possible preservation of n 

increases, where 

r 

[y,z]G0Q = gtj[yi,Zj]Q, y,z€X^, yi,...,yr€X, 

i,i=l 


and it will always be assumed that G G is Hermitian and non-singular. It 

is known m that the conditions for = [y[n-ll,y[n-l]]G^Q are that 

there exists a real diagonal s x s matrix D such that 


M : = 


DA + A^D - B^GB DU - B'^GV 
U^D-V^GB G-V^GV 


= 0 . 


(41) 


Note that A is assumed to remain real, but the other coefficient matrices may 
become complex—valued under a complex coordinate transformation T. Below, 
the method is assumed to be expressed in the canonical coordinates of Subsection 

ESI 
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Theorem 10 Let (^4, U, B, V) be a consistent method with real non-singular diagonal 
matrix D =diag(_B”ei). Then, the method is G-symplectic if and only if 

iV(C) = -D-^N^{C-^)D, C€C\Ao. (42) 


(Here, means evaluate the matrix function A'"” for the argument C~^-) 

Remark: Identity (I42II is a G-free characterization of G-symplecticity. In the 
linear multistep case [6], this is the same as (|38ll. 

Proof (only if) For € C \ (eig(F) U Z\o), fflll implies that 


0 = [ 7 , 




DA + A^D - B^GB DU - B^GV 
U^D-V^GB G-V^GV 


= [DA + DU{(I - Vy^B] + + B''(C"^7 - V'')" Vd] 

+(i - c-c~^)-B"(c"^ - f")“^g(c7 - vyy 

= 7)A1(C) +iV"(C“^)7). 


I 

{y-v)-y 


(43) 


(if) From (1391) and Lemma [5] it follows that there is a nonsingular T £ such 
that 


'A U' 


'D~^{B''V~''U'^ - A'^)D D~^B''V~'^T' 

B V 


rp-ly-Hj, 


From the (2,1) quadrant, B = T-^V~^U^D, and hence DUV~^ = B^T^. From 
the (1,2) and (2,2) quadrants, DUV~^ = B^T. Thus, 

0 = DUV~^ - DUV~^ = B^{T - T"). 

From the (2,2) quadrant, TV = V~^T, which implies T^V = Hence, 

4(r + r”)F = F””i(T + T"). 

Hence, \{T T”) may be substituted for T in (14411 : i.e. T may be assumed to be 

Hermitian. We now observe that (01) implies (01, with G = T. 


4.5 Methods that are both symmetric and G-symplectic 


Theorem 11 If a GLM satisfies two of the following conditions, it satisfies all three: 

(i) The method is symmetric; 

(ii) The method is G-symplectic; 

(ill) There exists a non-singular T £ 


'PDA PDU' 


' [PDAY B^T' 

TB TV 


[PDUY V^T 


(45) 


Condition (Hi) is equivalent to 


PDN{() = {PDN{())\ C€C\Zio. 


(46) 
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Proof The equivalence of (SSI) and (1461) follows from Lemma [21 On the other hand 
symmetry and G-symplecticity are respectively equivalent to the transfer func¬ 
tion indentities (1381) and (1421) . The following equivalences for the transfer function 
complete the proof: 

S6l)-hS21)^66l) + SHl)^SHl) + SI- 

The following closely connected result, the proof of which we omit, is useful 
in the construction of methods that are both symmetric and G-symplectic. The 
canonical coordinates of Subsection 13.51 are assumed. 

Theorem 12 Consider a consistent (L, P) -symmetric general linear method, where 
A is lower triangular and P is the reversing permutation matrix; i.e. {Pv)i = Vs+\-i, 
1 < i < s, for u € C®. Then, the method is G-symplectic if 

(i) non-zero real scalars hi, h 2 , ■ ■ ■ ,hr exist such that 

DUci = h^QB^'ei, I <i<r, (47) 

where hi is such that D = hi diag(B“ei). 

(ii) The diagonal part of A satisfies 


diag{aii,a 22 , • •., as-i,s-i, ass) = diag{ass, as-i,s-i, • • • , 022 , 011 )- (48) 

If the eigenvalues of V are distinct and diag(D) has no zero elements, then conditions 
(jiTl) and (1481) are also necessary for G-symplecticity. 


5 Symmetry and even order results 

5.1 Even order for the general linear method 

The method is of order p G N relative to the starting method Su if 


ShEhyo - MhShVo = Cp+i{yo)hP+^ + 0{hP+^), (49) 

where, for Tp+i the set of rooted trees of order p -|- 1, elementary differentials 
F(t)(po) £ X, symmetry coefficients a(t) € K and weight vectors !^'(t) G C*", 

iGTp+1 


The order of the method M/j is p G N, if p is the greatest integer such that there 
is an Sh relative to which 'Mh has order p. 

Following the work in Subsection 13.51 we assume that the method may be 
written in coordinates such that V, B and U take the form 


V = 


1 0 ^ 
0 V 



U=[1U]. 


(50) 


In particular, we note that B-i — E is non-singular. 

We assume that the method is of of order p relative to the starting method 
S/j. Written in the new basis, the principal component of is represented by 
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the B-series (see [E]). The remaining components are given by the vector of 
B-series, For some r;, representing the stage values, the stage equations and the 
update equations for the principal and non-principal components may be written 
in terms of B-series: 

rj{t) = A{r]D){t) + lC{t) + U^{t), (51) 

{ECm = b^{r,D){t) + (:{t), (52) 

{E^){t) = B{r^D){t) + Vm, (53) 

for all t such that |t| < p. Suppose that a second starting method Sh is similarly 
represented by B-series ( and 

Lemma 3 Suppose that the method M/j is of order p relative to and also of order 
p relative to S^, and that cg = cg. Then, 

C(i)=C(^)> |i|<P-l, 

fj{t)=ri{t), |t|<p-l, 

C(f)=C(^). \A<P- 

where rj is defined by (EH but for the starting method [CC]- 

Proof We first recall and extend some notation on trees. If |ti|,..., \tn\ > 1 then 

t = • • • tn] (54) 

denotes a rooted tree with order 

|t| = 1 + m -|- |tl I -|- \t2\ |tri| 

formed by joining the roots of m copies of r and each of the roots of ti {i = 
1,2,... ,n) to a new root. 

The valency of the root of t, will be written as 

w{t) = m + n. 

The binary product of trees will be used in the special case 

tr = ■ ■■tn], 

where t is given by (I54II . Note that wftr) = w{t) 4- 1. 

If p is the B-series representing stage values of a general linear method, then 
for this same t, the B-series for the stage derivatives are given by 

n 

{vD){t) = pir)""Y[piti), 

i=l 

where the powers and products on the right-hand side are componentwise. We will 
prove by induction on fc = 1, 2 ,... ,p — 1, 


II 

\t\ < k. 

(55) 

p{t) = p{f), 

|t| < k. 

(56) 

m = 

\t\ <k + l. 

(57) 
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Note that (I55II and (1561) are true when k = 0, and (i) follows from (1531) by substi¬ 
tuting the tree t = t to give 

with the same result for ^(r). Now assume the result for integers less than fe, 
and we prove (1551) for a specihc k € {1,2, 1}. For k = 1, this holds by 

assumption. For fc > 1, consider each tree t of order fc in a sequence in which w{t) 
is non-increasing. For t given by ([Ml). substitute tr into ([M]) to give the result 

n 

(m -h l)C(t) = J]^ r)(ti) - C{t, (), 

i=l 

where involves trees already considered for lower k and for trees with this 

same order which occurred earlier in the sequence. Obtain a similar result for 
and note that the terms on the right-hand side are identical in the two cases. The 
result (1561) follows from (1511) and the corresponding formula for To prove (15611 
for any tree of order k + 1, use (|53|) to obtain a formula for (/ — with the 

same result for (7 — 

Remarks: (i) The proof of LemmajSjserves as a constructive proof of the existence 
of S/j. Note that the order hP coefficient of C, is arbitrary. 

(ii) The assumption Ctt = CO can always be assumed because, if it were not true 
then S/j can be replaced by ShEeh for a suitable 0 £ K. 

Lemma 4 Suppose that the method is symmetric and of order p relative to 
such that CO “ 0- Then, M/j is also of order p relative to both L§_f, and the symmetric 
starting method ^{Sh+LS—h). The B-series for all3 starting methods agree up to order 
p, except possibly in the first component of the trees of order p. 

Proof Consider (1491) with h and yo replaced by —h and yi = Ef^yo- A left-multip¬ 
lication by 7/M“]j then yields 

{LMZlL){LS_h)E_hyi = LS.hVi + LV-^Cp+i{yi){-h)P+^ + 0{hP+^), 

where we note that the Frechet derivative of LJAz], is LV~^ + 0(h). Symmetry 
implies also, Cp+i{yi) = C'p+i(yo) + 0{h). Thus, 

Mh{LS_h)yo = {LS_h)E!,yo + LV-^Cp+i{yo){-h)P+^ + 0{hP+^), (58) 

and so M/j is of order p relative to L§_f,. Now, by LemmajSjthe B-series for 
and therefore also that for ^{Sh + LS-h), agrees with the B-series for S/j up to 
order p, except possibly in the hrst component of the trees of order p. The proof 
of LemmajSjshows that this is sufficient for -|- LS_h) to be a starting method 
relative to which M/j is of order p. 

Lemma 5 If (M/j, Sh) satisfy (1491) for some p £ N, then Sh may be chosen so that 
ShEhVO = ^hShVO + Kp+i{yo)eihP+^+0{hP+^), (59) 


where Kp+i{yo) := e{C'p+i(po) £ X- 
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Proof Replace §h satisfying (SHI), by Bh + Sh such that 

(7 - V)Sh = -{I - e^el)Cp+i{yo)hP+\ el5h = 0. 

Theorem 13 Suppose that M/j is a symmetric consistent method and that 1 is a 
simple eigenvalue ofV. Then, JAy is of even order p, and there is a symmetric starting 
method relative to which JAy is of order p. 

Proof Since JAy is consistent, it is of order p, for some p £ N. Lemma 2] ensures the 
existence of a symmetric starting method By relative to which M/j is of order p. 
Since LB_y = By, identities (14911 and (15811 are the same for this By. Equating the 
terms of order h^^^, we obtain 

Cp+i{yo) = {-l)PLV-^Cp+i{yo). 

By Lemma [5l Cp+i{yo) = Kp.^-i{yo)ei. As My is of order p, the term Kp^i{yo) is 
non-zero. From Subsection 13.51 LV~^ei = ei. Thus, 

Kp+i{yo) = elCp+i{yo) = {-lfelLV-^Cp+i{yo) = {-lfKp+i{yo). 

Hence, p is even. 

Lemma 6 Suppose that 1 is a simple eigenvalue ofV and that My is of order p relative 
to Biy. Then, there is a finishing method ly such that TyBy = I. If By = LB>_y, then 
Ty may be chosen to be symmetric; i.e. “Jy = T-yL. 

Proof Without loss of generality, we assume that By is represented as in Lemma 
|3l Since C(0) = 1, there exists an inverse B-series |121 11.12]. We observe that 
Hy := (~^el is a suitable hnishing method. Now suppose By is symmetric, and 
recall from Subsection 13.51 that Lei = ei and e\L = e{. Thus, LB-y = By implies 
that C is the same for $y and B_y. Hence, IP-yL = (^~^e\L = C~^e\ = “Jy. 

Theorem 14 Suppose that My is symmetric and of order p. Then, it is of order p 
relative to a symmetric starting method §y, with corresponding symmetric finishing 
method Hy. Furthermore, the error for initial data yo G X at x = nh is given by 

EnhVo - ^h’^h^hVo = h^Cp+iiyo, nh) + h^~'~^Cp+2{yo, nh) -, (60) 

where p is even and only even powers of h appear on the right-hand side of dMl). 

Proof The existence of suitable By and tpy is shown in Theorem [13] and Lemma 
|6| Let yo G X and a; £ K \ {0} be fixed. Given n £ Z \ {0}, define err(n) := 
Ex — 3 "2:/r!,^”/n^a:/rf Transforming n <—— n, and using the symmetry of My,By 
and !Ty, we obtain 

err(-n) = Ex - S'_,^/„MZ^/„B_,^/„ 

= Ex- T_,/„L(LM:1/„L)’^LS_,/„ 

Ex ^x/n^x/n^x/n err(n). 

Thus, err(n) is an even function of n. Hence, the expansion 

err(n) = (x/n)Pcp+i(yo, x) + (x/n)P~^^Cp+3(yo, *) H-, 

may only contain even powers of n. Putting h = x/n, we deduce that only even 
powers of h have non-zero coefficients in (16011 . 
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y[01 



Fig. 3: Commutative diagram for underlying one-step method 


5.2 The underlying one-step method 

Given a method 'Mhj the map <Pfi : X —> X is an underlying one-step method 
(UOSM) for M/i if there is a map : X —> X’’ such that 

^h^hVo = M/jS/jj/o, Vo £ X. (61) 

Relation (1611) may be represented by a commutative diagram as in Figure [31 

The concept of an underlying one-step method in the linear multistep case 
is due to Kirchgraber [16]. The existence of an underlying one-step method was 
extended to strictly stable general linear methods and made precise by Stoffer |25| . 
For the broader class of zero-stable methods, the existence and uniqueness of a 
formal B-series for and was shown in m- 

Because 'I'h = Ix + 0{h), : X —>• X is invertible and is also a 

UOSM for M,,: 


= M^(SA)yo, yo e X. (62) 

This freedom in and might be restricted in several ways. In m this is 
achieved by choosing a finishing method 'Jh : X’’ —> X in advance, and enforcing 
the finishing condition 

^hSh=Ix- (63) 

Here, we prefer to specify (, the B-series of the first component of S^- 

Below, we use the notation defined for LemmajSl and define B-series ip and [C, 
to represent and respectively. Equation (EH now implies the tree identities 

ri(t) = A(r/D)(t)+ l((t) + (64) 

{(l’Oii) = b'{vD){t) + C(t), (65) 

{<l,0{t) = B{vDm + Vm, ( 66 ) 

for a B-series rj representing the stage values. 

Theorem 15 Let JAh be a consistent zero-stable general linear method, such that the 
method may he written in the form (1501) with 1 a simple eigenvalue ofV. If C, is chosen 
such that (((0) = 1, then there exist unigue S/j and formally satisfying (16111 . 
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Proof For k = 0, (l6ll l65ll and (l66l) imply that ^(0) = 0, 77 ( 0 ) = 1 and 9s(0) = 1. 
For fc £ N, assume that (IMl 1651) and (l6^ hold for \t\ < k — 1. For |t| = k, r]{t) 
and are successively fixed by the following uniquely soluble rearrangements 
of (Ibh[M|l and (l65ll : 

m = {i- vr\B{nD){t) + e(t) - 

ri{t) = A{riD){t) + l({t) + U^{t), 

ip{t) = b'[iqD){t) + ((p(t) + C(t) - (pC)(i))- 

We observe that the terms on the right-hand side of the first and third equations 
depend only on the given value of (^(t) and on trees of order less than k. Once ^(t) 
is found, the second equation fixes r/(t). Induction on k now implies the existence 
of suitable rj and (p. Hence, there exist formal series for and Bf^ satisfying 
identity ([6T|). 

Remark: If ^ is chosen equal to the first component of the practical starting 
found in Lemma[3l then is a solution of (16511661) up to 0{hP). In that case, we 
deduce that the corresponding one-step method satisfies 


EhVO - ^hVo = Cp+i{yo)hP+^ + 0{hP+^). (67) 

Corollary 2 Let the assumptions of Theorem 1151 hold for a symmetric method Jily, 
and suppose that 

({t) = 0, |t| odd. (68) 

Let (Lfi and denote the corresponding underlying one-step and 

starting methods. Then, and are symmetric: 

= LB_h, <Th = (69) 


where denotes a formal inverse, and 

§(2«) ^ ^g(2g)^ g(2g-tl) ^ _^g(2qr+l)^ 5 = Q, 1, 2, . . . ( 70 ) 

Furthermore, dsa holds with p even. 


Proof Let {By, T’y) be as in the conclusion of Theorem 1151 Let h be replaced by 
—h in m and let yo be replaced by A left-multiplication by L then yields 

{LM_yL){LB_h)cl>zlyo = {LB_h)yo, 

(where all identities hold as formal B-series). Left-multiplication by M/j implies 
that 

{LB_h)L>-_lyo=Mh{LB_h)ye. 

Hence, {LB-y, ^Zy) s-lso satisfy (16111 . Now, by virtue of (1681) and Lei = ei, the 
B-series for the first component of LB-h is equal to C,, the first component of 
S/j. Thus, Theorem [TSl implies that {LB_h, <LZ\} = (S/j, ^/t), and we deduce (l6^ . 
Identities (1701) follow from a comparison of the coefhcients of h?'^ and in the 

expansions oi By and LB_h. 
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Substitute y-i = E_hyo for yo in (EZI) and left-multiply by Then, 

>l>.hyo = <P-hEhy-i = ^-h {'I'hy-i + Cp+i{E_yyo)hP+^ + 0{hP+^)'j 
= E_yyo + Cp+i{yo)hP+^ + 0{hP+^), 

where ‘P-h = Ix + 0{h). Under the transformation h <—> —h, we obtain 
Ehyo-'l>hyo = {-l)Pcp+iiyo)hP+^ +o{hP+^). 

A comparison with ([671) reveals that p is even. 


6 Examples of symmetric non—parasitic methods 

In this section we construct a number of symmetric methods, each of which is 
consistent and free of parasitism. Because we will consider only methods for which 
r = 2 and V = diag(l, —1), the parasitism there is only a single parasitism growth 
factor, equal to —e^BUe^., [2]. Parasitism growth rates are also discussed in m- 
For convenience, we select methods for which A is lower triangular, preferably 
with some zero elements on the diagonal. Many of the methods have r = 2 with 
V = diag(l, —1), and some are G-symplectic. For this choice of V, the two options 
L = diag(l, 1) and L = diag(l, —1) are possible and examples will be given for each 
of these. The terminology pqrs = 4123 indicates that there are r = 2 and s = 3, 
with order p = 4 and stage-order q = 1. Note that an irreducible method with 
rs = 22 can never be free of parasitism because for such a method, 622 = ±^21 and 
U 22 = ±'Ui 2 and hence the (2,2) element of BU equals 262 i'Ui 2 and this can only 
be zero if the method is reducible. Hence, we will start our examples with rs = 23. 


6.1 Starting and finishing methods 

We will present methods with r = 2 and L = diag(l, ±1). For ± = -|-, the principal 
input will be an even function and the second input will be an odd function. 
Suppose the B-series for these are defined by the coefficient vectors and ^ 2 , 
where 

ei(.) = = Ci(i) = ui) = = C2)(Y) = e2(j) = o, 

then it will be sufficient to also specify the required values of a; = [a;i, 0 : 2 , ® 3 , 2 : 4 ], 
where 

3:^1= C2(.), 3:2=^i(I), X 3 =^ 2 ^, 2:4 =^2(1). 

Note that the values of ^i{t) where |t| = 4 are irrelevant to the construction of 
appropriate starting values. Consider the two Runge-Kutta methods 


c 

A 

c 

A Pc-lb^l 

PAP - WP 


b'^’ 

V ■ 

-b'^P 


where P is the stage reversing permutation matrix. Note that the two Runge- 
Kutta methods are exact inverses. Hence if Ry is the mapping associated with 
{A,b^,c), then will be the mapping associated with {A,b^,c). 
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Impose on the {A, h', c) method the order conditions 

ClAl = XI, 

c = X2, 

Cb^ (? = *3, 

Cb^ Ac = Xi, 

where (7 is a constant at onr disposal. Based on and , we will use a starting 
method S^, defined by 

Vi^ = \{RhVQ + R-hVo), 

1 / 2 °' = i.C{Rhyo - R-hVo)- 

Similarly, we will use a finishing method F)^, defined by 

yn = \ [Rh{yY^ + U212”') + R-h{yY^ - ^1/2"'))• 

These proposed starting and finishing methods have the property that Ff^oSh = id 
and that they are consistent with the symmetry of the main method. 

Starting methods will be presented in the form of 

[C,iA,b\c),iA,V,c)]. (73) 


6.2 Methods with rs = 23 

Because we will insist on consistent, irredncible, parasitism-free methods, we will 
need to reject the case L = diag(l, 1). The reason for this is that symmetry wonld 
require 622 = 0, 623 = —621 and also M 22 = 0, M 32 = —mi 2 - Hence, the parasitism 
growth factor wonld he fi = — 2 fe 2 iwi 2 , and this would only be zero if either 621 = 0 
or M 12 = 0. However, in each of these cases, the method rednces to a Rnnge- 
Kutta method. However, methods exist with L = diag(l, —1) and the general case, 
assuming lower triangnlar A is given by 


On 

0 

0 

1 

Ml 

021 

022 

0 

1 

M2 

031 

032 

033 

1 

Ml 

bi 

62 

bi 

T" 

0 

./ 3 i 

/32 

Pi 

0 

-1 


snbject to 


261 + &2 = 1 , 

2/3imi -)- P 2 U 2 = 0, 

an + 033 = bi — ( 3 iui, 
0,21 = bi — /3iM2, 

2022 = &2 — h'02, 

031 = bi — Piui, 

032 = &2 — P2Ul- 
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By consistency, the methods in this family have order 1 and therefore, by Theorem 
El the order is also 2. For order 3, conditions associated with the trees of that 
order must be satisfied and, in this case, again using the even order result, the 
order must be 4. 

We present three examples of symmetric methods with rs = 23 and order 
4. None of these can be G-symplectic because this additional requirement would 
contradict the parasitism-free condition. 


First 4123 method 


The tableau for the method, which we will name 4123A, is 


A U 
B V 




0 0 

B 0 
2 1 

^4 


0 

-1 


To verify the order 4, we need to find a starting method, = S^yo such that 
the output after a single step of the method is = Shy{xo + h) + 0{h^). For this 
method a suitable choice of the starting values is given by 

= y{xo) + ^y"{xo), j/2°' = 

We note that = LS_h, as required for a symmetric starting method. We need 
to confirm that the result found by one step of the method is, to within 0{h^), 
equal to 

yW = y[xo + h) + ^y "(®0 + h) 

= y{xo) + hy{xo) + {xq) + 

2 ^ 2 ^' = (®o + h) - ^y"'i.xo + h) 

= ^y'{xo) + ^y"{xo) + '^y"'{xo) + ^y‘'‘^\xo). 

The B-series coefficients for and j/ 2 °^, corresponding to a tree t are denoted by 
and ^2 respectively, with the target values of the components of given by 
the components of E^. These are shown in Table [1] for the empty tree 0 and for 
the 8 trees of order up to 4. Also shown are the B-series coefficients for the three 
stages, denoted by ry and the stage derivatives {r]D)i, i = 1,2,3. Note that the 
table does not give values for %(t) where |t| = 4, because these are not needed in 
the evaluation of rjD up to order 4. Practical starting methods can be found in the 
form (17^ satisfying the order conditions for x = [ 5 , The solution is 


1 

4 

12, 0 

i 0 ^ 

4 24 

11 7 

12 12 ’ 24 

^ 0] 
1 1 
? ? 


1 1 

12 8 

1 1 

8 12 _ 


(74) 
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Table 1: Verification of the order of the method 412321 


aw 

a(t) 


Vl { t ) 

{viD){t) 

V 2 { t ) 

V 3 { t ) 
iV3 D)(t) 


(Eri 2 ){t) 


t V I v-vlY I 


10^ 0 
0 


1 0 J8 
0 10 

I 4 

2 4,8 

0 1 i 

II — 
Oil 


“TXT 

48 

0 

J_ 

f 

A. 

48 

1 


-fr 

1 

1 


0 0 0 
0 0 0 


^ 1 I 


0 0 Js 


1 _5_ 

8 96 48 


1 1 ^ _ _ 

' 1 ¥ 1 ¥ 

^ 2 2 16 32 


f : f f 


72 


96 


25 

144 


¥ ¥ f 


f 

96 


Here, as for the other methods in this section, one may choose the starting method 
to be explicit at the price of a more implicit finishing method, as the following 
alternative starting—finishing combinations indicate: 


4 

7’ 


0 

3 

I 

8 

7 

0 

3 
I 

1 

4 

0 

0 

25 

28 

0 

0 

0 ’ 

^|grH|oor-|oo 

1 1 

49 

648 

49 

Bis 

49 

648 

1 1 

O' 

0 

3 

4 


3 

65 

49 


49 

65 

3 


4 

324 

648 


648 

324 

4 J 


2VTE 

5 


Vl5 

12 

TvXs 

12 


SlyTS 

180 

169\/r5 

720 

31v^l5 

180 


4^/TE 

45 

aVTE 

45 ’ 

4v/l5 

45 


0 

0 

0 

%/l5 

16 

%/T5 

16 

0 


aivTE 

180 

4\/l5 
45 . 


Second 4123 method 


The following method, which we will denote as 41235, has the advantage of a zero 
on the diagonal. 


A U 
B V 


■i 0 0 

1-r 

Ml- 

O 

O 

1-s 

1 n 1 

1 “H 

ill 

3 3 3 

1 0 

1-2 1 

0 -1 


An analysis, similar to method 4123A, verifies order 4 with x = [Oj i^]- 

Although starting and hnishing methods similar to (Ell) do not exist, using two 
stage Runge-Kutta methods, they do exist with three stages. A possible triple is: 


0 

0 

0 

0 

1 

1 

23 

0 

V4 

24 

3 

i 

i 

2 

4 

24 

24 

3 




1 


24 

24 

12 


3 

4 

3 

4 

0 

0 

1 

1 

11 

0 


12 

12 

0 

1 

1 

1 

12 

24 

24 



1 

i 


12 

24 

24 . 


-12, 
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A 4223 method 

The following method, named 4223 A, is fonnd to have stage order 2 , 


A U 
B V 


■i 0 0 

i-i' 

0 i 0 

1 1 

13 1 
4 4 8 

1-^ 

12 1 
6 3 6 

1 0 

111 

.666 

0 -1 


Suitable starting values are 


2 / 1 °' = v(xo), 


[ 0 ] ^ f / \ m / \ 

y2 =jy (®o) - -^y (a:o), 


corresponding to a: = [- 1 , 0 ,—No hnishing method is required other than 
2/n = 2/1"' and the starting method can be defined by 1 / 2 °' = \{Rhyo — R-hyo) where 
Rh is the Runge-Kutta method with tableau 


0 


1 

1 

4 

4 

1 

4 

0 ^ 


7 11 


12 6 6 


A special 4123 method 

The method to be named 4123 C is defined by 


A U' 
B V 


"000 
X A 0 

12 12 ^ 

111 
12 6 12 

1 1 

1 -1 

1 1 


1 1 f 

3 3 3 

111 
-4 2 4 

1 0 

0 -1 


This method is interesting because, although it is symmetric, the diagonal of A is 
not symmetric. 

Using X = [5, —a starting-hnishing triple is found: 


12 , 


0 

0 

0 

1 

5 

1 

_ 4 _ 

24 

24 • 


b 

i 


24 

6 


5 

5 

0 ■ 

24 

24 

1 

1 

5 

24 


24 


1 

5 


6 

24 . 
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6.3 Methods with rs = 24 
First method with L = diag(l, —1) 

We now search for symmetric methods of the form 






1 

Ul 





1 

U2 



A 


1 

U2 





1 

Ul 

61 

^2 

b 2 

bi 

1 

0 

La 

P2 

P 2 

Pi 

0 

-1 


with P\Ui + P 2 U 2 = 0 (to eliminate parasitism) and order 4. We give an example 
which will be named 4124A: 



This method has the same symmetry, defined by L = diag(l, —1), as in Subsection 
16.21 and it is possible to use similar starting and finishing methods, An analysis, 
which will not be included, leads to a starting-finishing pair defined from x = 
[0, — — j^]. The triple defining the pair is 


0 

0 0 i 

1 0 ] 

18, i 

6 12 ’ 0 

1 1 

6 6 


i i 

i 1 

. 

6 6 

6 6 J 


G-symplectic method with L = diag(l, —1) 


By imposing the requirements of Theorem 1121 a G-symplectic symmetric method 
can be constructed with G = diag(l,—^), D = diag(—g, |, |, —^) and order 4. 
This method, denoted by 4124B, has the tableau 



is given by 


27^, 


0 

0 

0 

x/38 

5V^ 

%/38 

24 

152 

114 


5V^ 

V38 


152 

38 


2x/38 

2x/38 

0 

57 

57 


Vss 

a/38 

5a/38 

152 

38 

152 


a/38 

V38 


38 

152 . 
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Fourth order symmetric methods with L = diag(l, 1) 

We will derive symmetric parasitism-free methods with L = I and an = 0, based 
on the assumptions 


6^1 = 1, 


lT 2 1 

& C = 2 , 

Ac = g. 


(75) 

(76) 

(77) 


From m and dZSl), it is found that 

^2 = T7^ — 7Z -r, bi = - - 62. 

12c2(1-C2) 2 

Without loss of generality, because we can use a diagonal scaling transformation, 
assume mi = 1 and, to eliminate parasitism, it follows that /?i = —« 2 /? 2 - We will 
impose the condition an = an = 0, implying that /3i = 61 . From A + PAP = 
UV~^B, and the requirement that A is lower triangular we find that 


A = 


0 0 0 O' 

61 ( 1 -U 2 ) j+x 0 0 

61(1 +U 2 ) 62 -&i i-x 0 ’ 

2b^ b2-^b2 + ^0_ 


where x is arbitrary. The value of U 2 is determined by the requirement that Al = c 
and this gives 


_ -12ci + 21 ci- 9 c 2 + l 12 c 2 (c 2 - l)x 
^ 6 c| — 6 c 2 + 1 6 c| — 6 c 2 + 1 

For (1771) to be satisfied, a complicated condition is obtained. This is satisfied for 
any value of x if and only if C 2 = ^ and this is the value that will be selected. We 
present the matrices defining the method in three cases a; = — ^, a: = 0 and x = 
We denote the corresponding methods as 4124C, 4124D and 4124E: 




0 

0 

0 

0 

1 1 



i 0 0 0 

1 -2 

A U' 


_i i i 0 

6 6 2 ^ 

1 2 

B V 


1 5 1 r. 

3 12 1 

1 -1 



i i i 1 

6 3 3 6 

1 0 



1111 

0 1 



L 6 12 12 6 
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Each of these three methods has order 4 for identical conditions on the starting 
method. These are defined by 


t 

0 

• 

I 

V I 

V 


Y 

1 

Ci(i) 

1 

0 

0 

0 0 

0 

0 

0 

0 

?2(t) 

0 

0 - 

1 

12 

0 0 

1 

24 

1 

36 

1 

48 

1 

144 


Because corresponds to the identity mapping, the finishing method can be 
chosen as 

A practical starting method is available in the form = yo and = 
RhVo — yo, where is defined by the Runge—Kutta tableau 


0 




1 

1 



2 

2 



1 

5 

1 


2 

6 

3 


1 

4 

5 

1 

3 

6 

2 


I 

0 

1 1 


4 

3 12 


7 Simulations 

We compare the long-time numerical behaviour of several symmetric general lin¬ 
ear methods with that of two symmetric Runge-Kutta methods. One of these 
RK methods is symplectic, and two of the GLMs are G-symplectic. The four 
low-dimensional Hamiltonian test problems we consider have one or more of the 
following properties: absence of symmetry, non-separability, chaotic behaviour, or 
large time derivatives. We compare the efficiency of the methods, as well as their 
ability to conserve invariants over long times. 


7.1 The problems 

Henon-Heiles The equations of motion are defined by the separable Hamiltonian 
H{p,q) = ^{pi +pl) + +qi) + qiql - \q 2 - 
The initial conditions are taken [10] so that H = j ■. 

[pi,P2,?i, 92 ]''= 0 . 2 , 0 , 0 .dj . 

The solution is chaotic. In the experiments, the time-step h = 0.25, and the final 
time T = 10®. 
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Double pendulum The equations of motion are defined by the non—separable Hamil¬ 
tonian 


H {p, q) 


pI + 2 p 2 - 2 pip 2 cos(gi - q 2 ) 
2(1 -I- sin^(gi - 52 )) 


cos((j 2 ) — 2cos{gi). 


For y ■= [p'^,q'^Y, let f{y) := J ^WH{y), where 


J : = 


0 I 
-I 0 


0 0 10 

0 0 0 1 

-1 0 0 0 

0-100 


For R = diag(—1, —1,1,1) or diag(l, 1, —1, —1), the system is p-reversible [2l]; i.e., 
f{Ry) = -Rf{y), for y g 


The initial conditions are taken to be 

[pi,P 2 ,</i,<? 2 ]''= [o,0,3.14,-3.lj . 

The solution is chaotic with large time-derivatives. Here, h = 0.01 and T = 10^. 


Kepler problem This describes the motion of a planet revolving around the sun, 
which is considered to be fixed at the origin. The equations of motion are defined 
by the separable Hamiltonian, 

H{p,q) = ^{pI +pI) - , \ , 

Vgf + <ii 

where q = [ 91 , 52 ]^ are the generalized position coordinates of the body and p = 
[pi,P 2 Y are the generalized momenta. For y := [p'^,q'^Yj fof fiv) ■= J~^^H{y). 
Then, the system is multiply p-reversible for 

R = diag(-l, -1,1,1), diag(l, 1, -1, -1), diag{l, -1, -1,1), or diag(-l, 1,1, -1). 
The initial conditions are taken to be 

[pi,P 2 ,< 7 i,<? 2 ]'^ = [ 0 , 1 - e,oj , 

for e = 0.6. The solution is a closed orbit with moderately large time derivatives. 
The angular momentum error is plotted in addition to the Hamiltonian error. 
Here, h = 0.01 and T = 10^. 


Transformed Lotka-Volterra The equations of motion are defined by the separable 
Hamiltonian 

H{p,q) =p-exp(p)-|-2g-exp(g). 

This system lacks any obvious symmetry. The initial conditions are taken to be 

[p,qY = [ln2,ln3j . 

The solution is a non-symmetric orbit in the positive quadrant p, q > 0. Here, 
h = 0.1 and T = 10^. 
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Table 2: Timings for specific methods with various problems 




H-H 

DP 

K 

TLV 

4223 

Stopwatch 

CPU 

5.7160 X 10^ s 
5.7846 X 10® s 

1.5971 X 10^ s 
1.6153 X 10^ s 

1.0489 X 10^ s 
1.0629 X 10^ s 

8.8033 s 
8.7517 s 

4124B 

Stopwatch 

CPU 

7.4698 X 10^ s 
7.5623 X 10® s 

2.1707 X 10^ s 
2.1830 X 10^ s 

1.4358 X 10^ s 
1.4512 X 10® s 

11.7559 s 
11.8249 s 

4124D 

Stopwatch 

CPU 

4.9182 X 10^ s 
4.9746 X 10® s 

1.4974 X 10^ s 
1.4908 X 10^ s 

0.9398 X 10^ s 
0.9512 X 10® s 

7.9902 s 
7.9405 s 

5-DIRK 

Stopwatch 

CPU 

9.3599 X 10^ s 
9.9014 X 10® s 

2.5711 X 10^ s 
2.7857 X 10^ s 

1.5986 X 10^ s 
1.6254 X 10® s 

13.1807 s 
13.8061 s 


7.2 Methods used in the simulations 

The following methods are competitively compared in the initial simulations: 

— Method 4223 from Subsection 6.2: this is symmetric. 

— Method 4124B from Subsection 6.3: this is symmeric and G-symplectic. 

— Method 4124D from Subsection 6.3: this is symmetric. 

— The DIRK 4115 method: a 5-step Suzuki composition of the implicit midpoint 
2111 method, see [ISl Chapter II]: this is symmetric and symplectic. (This is 
more efficient than the familiar 3-step 4113 DIRK composition method due to 
far smaller error constants.) 

Simulations with two other methods serve to interpret the initial results: 

— The 4113 Lobatto IIIB method, [101 Chapter XI]: this symmetric, but not 
symplectic. 

— Method 4124P from [2]: this is symmetric and G-symplectic. 


7.3 Numerical simulations 

As with long-time Runge-Kutta experiments, we use compensated summation 
and a tight error tolerance for implicit iterations in an attempt to reduce the 
effects of rounding error. In order to reduce potential parasitic effects, we also use 
an accurate starting method for the multivalue experiments. 

Timings In Table [2] details of the CPU and stopwatch times for each of the exper¬ 
iments are summarised. 


7.4 Interpretation of the simulations 

Numerical errors in computing the invariants proceed from several potential sources: 

(a) The underlying one-step method does not possess the geometric properties re¬ 
quired for the problem. 

(b) Small periodic deviations occur, for example, when a conjugate symplectic 
UOSM approximately conserves a modified Hamiltonian Tf/j, and deviates 
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(a) 4223 



(b) 4124B 



(c) 4124D 



(d) 5-jump DIRK 
Fig. 4: Henon-Heiles problem 
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(a) 4223 


10® 10' 10^ 

0 - 

-3.11e-6 - 
-6.21e-6 - 
^ -9.32e-6 - 
-1.24e-5 - 
-1.55e-5 - 
-1.86e-5- 


steps 

10^ 


10'^ 


io’ 

time 


(b) 4124B 






(d) 5-jump DIRK 

Fig. 5: Double pendulum problem 
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steps 

10^ 



(d) 5-jump DIRK 

Fig. 7: Kepler problem: Hamiltonian 
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(b) Double Pendulum 


steps 

10“ 10’ 10^ 10^ 10^ 



(c) Transformed Lotka Volterra 

Fig. 9: Lobatto IIIB Experiments 


from the true Hamiltonian by a small, roughly periodic quantity. 

(c) Parasitism. 

Classically, we think of the effects of (a) and (c) as being clear-cut. However, 
the lack of symplecticity in high-order symmetric methods may take a very long 
time to manifest itself, see e.g. the behaviour of Lobatto HIA in [7]. This is also true 
of the effects of higher-order parasitism. In order to distinguish computationally 
the effects due to these two possible causes for the purely symmetric methods 
4223 and 4124D, we have also presented results for the 4113 Lobatto IIIB method, 
which has similar properties to the UOSMs of 4223 and 4124D. Finally, we have 
also shown results for the symmetric G-symplectic 4124P method applied to the 
TLV problem, as an improvement on those of 4124B. 
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steps 



(a) Kepler: Angular Momentum 



(b) Kepler: Hamiltonian 



Fig. 11: 4124P: Transformed Lotka-Volterra 


Henon-Heiles All methods exhibit broadly similar conservation behaviour. In the 
absence of parasitism, it is unsurprising that the results for the G-symplectic 
method 4124B should resemble those of the Suzuki 4115 DIRK. Also, the be¬ 
haviour of the purely symmetric 4223 and 4124D methods may be explained in 
terms of their UOSMs, which are closely related to the 4113 Lobatto IIIB method. 
Following the explanation of [7] for symmetric Runge-Kutta methods, the fact 
that H{p, q) is a cubic polynomial implies that the bushy trees in the numerical 
modified Hamiltonian vanish for order greater than 4. This permits the existence 
of an exact modified Hamiltonian for the UOSM of a symmetric non-symplectic 
method of order 4. Thus, even for chaotic solutions, one can expect conservation 
of a modified Hamiltonian, in the absence of parasitism. 
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Double pendulum All methods exhibit broadly similar conservation behaviour. The 
system is p-reversible, but as the behaviour is chaotic, no analog of the symmet¬ 
ric conservation result, [ini Theorem XI.3.1], would seem to hold in this case. 
Comparing the graphs for 4124D and 4113 Lobatto IIIB, we see broadly similar 
behaviour. We would therefore attribute any minor deviations in the Hamiltonian 
as due to properties of the UOSM, rather than to higher-order parasitism. 

Kepler The quadratic angular momentum is exactly conserved by the symplectic 
Suzuki 4115 DIRK, apart from random round-off errors. Otherwise, all methods 
exhibit similar conservation behaviour. Again, in the absence of parasitism, this is 
what one would expect for the G-symplectic 4124B method. The conservation be¬ 
haviour for 4223 and 4124D follows that of 4113 Lobatto IIIB. In this case, Kepler 
is both integrable and reversible. Although the exact hypotheses of [B Theorem 
XI.3.1] are not satished here, the situation is sufhciently similar to conjecture that 
symmetric UOSMs conserve invariants to 0(hP) uniformly in time, in the absence 
of parasitism. 

Transformed Lotka-Volterra This is a Hamiltonian problem without symmetry. In 
the initial simulations, only the Suzuki 4115 DIRK exhibits satisfactory approxi¬ 
mate conservation of the Hamiltonian. The lack of symmetry in the problem and 
the lack of symplecticity in the UOSMs for 4223 and 4124D methods explains the 
poor results in those cases. Although 4124B roughly conserves the Hamiltonian, 
there is a hint of parasitism at the end of the computation. The results for the 
G-symplectic method 4I24P show that good conservation is possible for general 
linear methods. 


7.5 Conclusions 

All methods performed similarly on the first three problems: Henon-Heiles, Dou¬ 
ble Pendulum and Kepler, except that angular momentum was exactly conserved 
only by the exactly symplectic Runge-Kutta method. Although the errors for the 
Suzuki 4115 DIRK were about 4 times smaller than those of 4124D for the hxed 
time-steps used, the timings indicate that the latter method is slightly more ef¬ 
ficient. Since 4124D only has 2 implicit stages, one would expect this efficiency 
advantage to increase for larger problems. 

Although parasitism did not develop for these problems, despite chaotic be¬ 
haviour, large derivatives and long time—intervals, further theoretical work and 
computational tests would be needed before general linear methods could be ap¬ 
plied to other problems with complete confidence. In the absence of parasitism, 
it appears that symmetric general linear methods behave in the same way as 
symmetric Runge-Kutta methods, whilst G-symplectic GLMs behave similarly to 
symplectic RKMs, with the exception that quadratic quantities are not conserved 
exactly. In particular, symmetric GLMs are not suitable for non-symmetric Hamil¬ 
tonian systems, such as the transformed Lotka-Volterra problem. 
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